
!------------------------------------------------------------------------!
!  The Community Multiscale Air Quality (CMAQ) system software is in     !
!  continuous development by various groups and is based on information  !
!  from these groups: Federal Government employees, contractors working  !
!  within a United States Government contract, and non-Federal sources   !
!  including research institutions.  These groups give the Government    !
!  permission to use, prepare derivative works of, and distribute copies !
!  of their work in the CMAQ system to the public and to permit others   !
!  to do so.  The United States Environmental Protection Agency          !
!  therefore grants similar permission to use the CMAQ system software,  !
!  but users are requested to provide copies of derivative works or      !
!  products designed to operate in the CMAQ system to the United States  !
!  Government without restrictions as to use by others.  Software        !
!  that is used with the CMAQ system but distributed under the GNU       !
!  General Public License or the GNU Lesser General Public License is    !
!  subject to their copyright restrictions.                              !
!------------------------------------------------------------------------!

       SUBROUTINE HRG4( DTC )

C**********************************************************************
C
C  FUNCTION:  To solve for the concentration of NO3 and N2O5
C
C  PRECONDITIONS: For the CRACMM2 mechanism
C
C  KEY SUBROUTINES/FUNCTIONS CALLED: None
C
C  REVISION HISTORY: Created by EBI solver program, Mar 15, 2024
C
C   18 Jul 14 B.Hutzell: revised to use real(8) variables
C**********************************************************************
      USE HRDATA

      IMPLICIT NONE


C..INCLUDES: NONE


C..ARGUMENTS:
      REAL( 8 ), INTENT( IN ) :: DTC             ! Time step


C..PARAMETERS: NONE


C..EXTERNAL FUNCTIONS: NONE


C..SAVED LOCAL VARIABLES:
!     CHARACTER( 16 ), SAVE   ::  PNAME = 'HRG4'    ! Program name


C..SCRATCH LOCAL VARIABLES:
      REAL( 8 ) ::   A, B, C, Q   ! Quadratic equation terms
      REAL( 8 ) ::   CMN          ! Temp scalar
      REAL( 8 ) ::   L15          ! Loss of NO3
      REAL( 8 ) ::   L16          ! Loss of N2O5
      REAL( 8 ) ::   P15          ! Production of NO3
      REAL( 8 ) ::   K15_15       ! Kno3+no3 * delta t
      REAL( 8 ) ::   R15_16       ! Kn2o5-->no3 * delta t
      REAL( 8 ) ::   R16_15       ! Kno3+no2-->n2o5[NO2] * delta t


C**********************************************************************


c..Production of NO3 (except from N2O5 )
      P15 =    2.0000D-01 * RKI(     9 ) * YC ( HNO4        )                           ! HNO4=0.2000D+00*NO3+...
     &    +                 RKI(    38 ) * YC ( PAN         )                           ! PAN=NO3+MO2
     &    +                 RKI(    44 ) * YCP( O3          ) * YCP( NO2         )      ! O3+NO2=NO3
     &    +                 RKI(    62 ) * YCP( NO2         ) * YCP( O3P         )      ! NO2+O3P=NO3
     &    +                 RKI(    64 ) * YC ( HNO3        ) * YCP( HO          )      ! HNO3+HO=NO3
     &    +                 RKI(   138 ) * YC ( PAN         ) * YCP( HO          )      ! PAN+HO=NO3+XO2+HCHO
     &    +                 RKI(   139 ) * YC ( PPN         ) * YCP( HO          )      ! PPN+HO=NO3+XO2+HCHO
     &    +                 RKI(   143 ) * YC ( HONIT       ) * YCP( HO          )      ! HONIT+HO=NO3+HKET

c..Loss frequency of NO3 ( except NO3 + NO3 if present )
      L15 =                 RKI(     5 )                        ! NO3=NO
     &    +                 RKI(     6 )                        ! NO3=O3P+NO2
     &    +                 RKI(    65 ) * YCP( HO          )   ! NO3+HO=HO2+NO2
     &    +                 RKI(    66 ) * YCP( HO2         )   ! NO3+HO2=0.7000D+...
     &    +                 RKI(    67 ) * YCP( NO          )   ! NO3+NO=0.2000D+...
     &    +                 RKI(    68 ) * YCP( NO2         )   ! NO3+NO2=NO+NO2
     &    +                 RKI(    70 ) * YCP( NO2         )   ! NO3+NO2=N2O5
     &    +                 RKI(    91 ) * YC ( ISO         )   ! NO3+ISO=0.4000D+...
     &    +                 RKI(   158 ) * YC ( ETE         )   ! NO3+ETE=0.8000D+...
     &    +                 RKI(   159 ) * YC ( OLT         )   ! NO3+OLT=0.4300D+...
     &    +                 RKI(   160 ) * YC ( OLI         )   ! NO3+OLI=0.1100D+...
     &    +                 RKI(   161 ) * YC ( API         )   ! NO3+API=0.9750D+...
     &    +                 RKI(   162 ) * YC ( LIM         )   ! NO3+LIM=0.9450D+...
     &    +                 RKI(   163 ) * YC ( TRPN        )   ! NO3+TRPN=0.3300D+...
     &    +                 RKI(   164 ) * YC ( HCHO        )   ! NO3+HCHO=HO2+CO+HNO3
     &    +                 RKI(   165 ) * YC ( ACD         )   ! NO3+ACD=ACO3+HNO3
     &    +                 RKI(   166 ) * YC ( ALD         )   ! NO3+ALD=RCO3+HNO3
     &    +                 RKI(   167 ) * YC ( MACR        )   ! NO3+MACR=0.6800D+...
     &    +                 RKI(   168 ) * YC ( UALD        )   ! NO3+UALD=HO2+XO2+...
     &    +                 RKI(   169 ) * YC ( GLY         )   ! NO3+GLY=HO2+...
     &    +                 RKI(   170 ) * YC ( MGLY        )   ! NO3+MGLY=ACO3+CO+...
     &    +                 RKI(   171 ) * YC ( PHEN        )   ! NO3+PHEN=0.1520D+...
     &    +                 RKI(   172 ) * YC ( CSL         )   ! NO3+CSL=0.2000D+...
     &    +                 RKI(   173 ) * YC ( MCT         )   ! NO3+MCT=MCTO+HNO3
     &    +                 RKI(   174 ) * YC ( MPAN        )   ! NO3+MPAN=MACP+NO2
     &    +                 RKI(   339 ) * YC ( MO2         )   ! NO3+MO2=HO2+HCHO+NO2
     &    +                 RKI(   340 ) * YC ( ETHP        )   ! NO3+ETHP=HO2+NO2+ACD
     &    +                 RKI(   341 ) * YC ( HC3P        )   ! NO3+HC3P=0.2540D+...
     &    +                 RKI(   342 ) * YC ( HC5P        )   ! NO3+HC5P=0.4880D+...
     &    +                 RKI(   343 ) * YC ( ETEP        )   ! NO3+ETEP=HO2+NO2+...
     &    +                 RKI(   344 ) * YC ( OLTP        )   ! NO3+OLTP=0.4700D+...
     &    +                 RKI(   345 ) * YC ( OLIP        )   ! NO3+OLIP=0.8600D+...
     &    +                 RKI(   346 ) * YC ( BENP        )   ! NO3+BENP=NO2+HO2+...
     &    +                 RKI(   347 ) * YC ( TOLP        )   ! NO3+TOLP=NO2+...
     &    +                 RKI(   348 ) * YC ( XYLP        )   ! NO3+XYLP=NO2+...
     &    +                 RKI(   349 ) * YC ( EBZP        )   ! NO3+EBZP=NO2+...
     &    +                 RKI(   350 ) * YC ( ISOP        )   ! NO3+ISOP=HO2+NO2+...
     &    +                 RKI(   351 ) * YC ( APIP1       )   ! NO3+APIP1=NO2+...
     &    +                 RKI(   352 ) * YC ( LIMP1       )   ! NO3+LIMP1=HO2+...
     &    +                 RKI(   353 ) * YC ( APINP1      )   ! NO3+APINP1=...
     &    +                 RKI(   354 ) * YC ( LIMNP1      )   ! NO3+LIMNP1=...
     &    +                 RKI(   355 ) * YC ( ACO3        )   ! NO3+ACO3=MO2+NO2
     &    +                 RKI(   356 ) * YC ( RCO3        )   ! NO3+RCO3=ETHP+NO2
     &    +                 RKI(   357 ) * YC ( ACTP        )   ! NO3+ACTP=ACO3+...
     &    +                 RKI(   358 ) * YC ( MEKP        )   ! NO3+MEKP=0.6700D+...
     &    +                 RKI(   359 ) * YC ( KETP        )   ! NO3+KETP=HO2+NO2+...
     &    +                 RKI(   360 ) * YC ( MACP        )   ! NO3+MACP=HCHO+...
     &    +                 RKI(   361 ) * YC ( MCP         )   ! NO3+MCP=NO2+HO2+...
     &    +                 RKI(   362 ) * YC ( MVKP        )   ! NO3+MVKP=0.3000D+...
     &    +                 RKI(   363 ) * YC ( UALP        )   ! NO3+UALP=HO2+NO2+...
     &    +                 RKI(   364 ) * YC ( BALP        )   ! NO3+BALP=BAL1+NO2
     &    +                 RKI(   365 ) * YC ( BAL1        )   ! NO3+BAL1=BAL2+NO2
     &    +                 RKI(   366 ) * YC ( ADDC        )   ! NO3+ADDC=HO2+NO2+...
     &    +                 RKI(   367 ) * YC ( MCTP        )   ! NO3+MCTP=NO2+MCTO
     &    +                 RKI(   368 ) * YC ( ORAP        )   ! NO3+ORAP=NO2+GLY+HO2
     &    +                 RKI(   369 ) * YC ( OLNN        )   ! NO3+OLNN=HO2+NO2+...
     &    +                 RKI(   370 ) * YC ( OLND        )   ! NO3+OLND=0.2000D+...
     &    +                 RKI(   371 ) * YC ( ADCN        )   ! NO3+ADCN=0.2000D+...
     &    +                 RKI(   375 ) * YC ( XO2         )   ! NO3+XO2=NO2
     &    +                 RKI(   395 ) * YC ( ACRO        )   ! NO3+ACRO=0.6800D+...
     &    +                 RKI(   399 ) * YC ( BDE13P      )   ! NO3+BDE13P=HO2+...
     &    +                 RKI(   404 ) * YC ( BDE13       )   ! NO3+BDE13=...
     &    +                 RKI(   410 ) * YC ( FURAN       )   ! NO3+FURAN=NO2+...
     &    +                 RKI(   412 ) * YC ( SESQ        )   ! NO3+SESQ=SESQNRO2
     &    +                 RKI(   415 ) * YC ( SESQNRO2    )   ! NO3+SESQNRO2=...
     &    +                 RKI(   419 ) * YC ( SESQRO2     )   ! NO3+SESQRO2=...
     &    +                 RKI(   426 )                        ! NO3=HNO3
     &    +                 RKI(   447 ) * YC ( VROCP6ALKP  )   ! NO3+VROCP6ALKP=...
     &    +                 RKI(   448 ) * YC ( VROCP5ALKP  )   ! NO3+VROCP5ALKP=...
     &    +                 RKI(   449 ) * YC ( VROCP4ALKP  )   ! NO3+VROCP4ALKP=...
     &    +                 RKI(   450 ) * YC ( VROCP3ALKP  )   ! NO3+VROCP3ALKP=...
     &    +                 RKI(   451 ) * YC ( VROCP2ALKP  )   ! NO3+VROCP2ALKP=...
     &    +                 RKI(   452 ) * YC ( VROCP1ALKP  )   ! NO3+VROCP1ALKP=...
     &    +                 RKI(   453 ) * YC ( HC10P       )   ! NO3+HC10P=HC10P2+NO2
     &    +                 RKI(   475 ) * YC ( VROCP6ALKP2 )   ! NO3+VROCP6ALKP2=...
     &    +                 RKI(   476 ) * YC ( VROCP5ALKP2 )   ! NO3+VROCP5ALKP2=...
     &    +                 RKI(   477 ) * YC ( VROCP4ALKP2 )   ! NO3+VROCP4ALKP2=...
     &    +                 RKI(   478 ) * YC ( VROCP3ALKP2 )   ! NO3+VROCP3ALKP2=...
     &    +                 RKI(   479 ) * YC ( VROCP2ALKP2 )   ! NO3+VROCP2ALKP2=...
     &    +                 RKI(   480 ) * YC ( VROCP1ALKP2 )   ! NO3+VROCP1ALKP2=...
     &    +                 RKI(   481 ) * YC ( HC10P2      )   ! NO3+HC10P2=NO2+...
     &    +                 RKI(   492 ) * YC ( VROCP6AROP  )   ! NO3+VROCP6AROP=...
     &    +                 RKI(   498 ) * YC ( VROCP5AROP  )   ! NO3+VROCP5AROP=...
     &    +                 RKI(   504 ) * YC ( NAPHP       )   ! NO3+NAPHP=NO2+...
     &    +                 RKI(   529 ) * YC ( STYP        )   ! NO3+STYP=NO2+HO2+...

c..Loss frequency of N2O5
      L16 =                 RKI(    71 )                        ! N2O5=NO2+NO3
     &    +                 RKI(    72 )                        ! N2O5=0.2000D+01*HNO3
     &    +                 RKI(   423 )                        ! N2O5=0.2000D+01*HNO3

c..K15_15, R15_16, and R16_15 terms
      K15_15  = RKI(    69 ) * DTC

      R15_16  = ( RKI(    71 ) ) * DTC 


      R16_15  = RKI(    70 ) * YCP( NO2 ) * DTC

c..Solution of quadratic equation to get NO3 & N2O5
      CMN = 1.0D0 + L16 * DTC
      A = 2.0D0 * K15_15 * CMN
      B = CMN * ( 1.0D0 + L15 * DTC ) - R15_16 * R16_15
      C = CMN * ( YC0( NO3 ) + P15 * DTC ) +  R15_16 * YC0( N2O5 )

      Q = -0.5D0 * ( B + SIGN( 1.0D0, B ) * SQRT( B * B + 4.0D0 * A * C ) )
      YCP( NO3 ) = MAX( Q / A , -C / Q  )
      YCP( N2O5 ) = ( YC0( N2O5 ) + R16_15 * YCP( NO3 ) ) / CMN

      RETURN

      END
